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ABSTRACT 

We present the angular power spectrum derived from the first-year Wilkinson Microwave 
Anisotropy Probe (WMAP) sky maps. We study a variety of power spectrum estimation meth- 
ods and data combinations and demonstrate that the results are robust. The data are modestly 
contaminated by diffuse Galactic foreground emission, but we show that a simple Galactic tem- 
plate model is sufficient to remove the signal. Point sources produce a modest contamination in 
the low frequency data. After masking ^700 known bright sources from the maps, we estimate 
residual sources contribute ^3500 /iK 2 at 41 GHz, and ^130 /iK 2 at 94 GHz, to the power spec- 
trum [1(1 + l)Ci/2n] at I — 1000. Systematic errors are negligible compared to the (modest) level 
of foreground emission. Our best estimate of the power spectrum is derived from 28 cross-power 
spectra of statistically independent channels. The final spectrum is essentially independent of 
the noise properties of an individual radiometer. The resulting spectrum provides a definitive 
measurement of the CMB power spectrum, with uncertainties limited by cosmic variance, up to 
I ~ 350. The spectrum clearly exhibits a first acoustic peak at I = 220 and a second acoustic 
peak at I ~ 540 (Page et al. 2003b), and it provides strong support for adiabatic initial condi- 
tions (Spergel et al. 2003). Kogut et al. (2003) analyze the Cj E power spectrum, and present 
evidence for a relatively high optical depth, and an early period of cosmic reionization. Among 
other things, this implies that the temperature power spectrum has been suppressed by ~30% 
on degree angular scales, due to secondary scattering. 
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1. 



INTRODUCTION 



The Wilkinson Microwave Anisotropy Probe (WMAP) mission was designed to measure the CMB 
anisotropy with unprecedented precision and accuracy on angular scales from the full sky to several arc 
minutes by producing maps at five frequencies from 23 to 94 GHz. The WMAP satellite mission (Bennett 
ct al. 2003a) employs a matched pair of 1.4m telescopes (Page et al. 2003c) observing two areas on the sky 
separated by ~141°. A differential radiometer (Jarosik et al. 2003a) with a total of 10 feeds for each of the 
two sets of optics (Barnes et al. 2003; Page et al. 2003c) measures the difference in sky brightness between 
the two sky pixels. The satellite is deployed at the Earth-Sun Lagrange point, L 2} and observes the sky 
with a compound spin and precession that covers the full sky every six months. The differential data are 
processed on the ground to produce full sky maps of the CMB anisotropy (Hinshaw et al. 2003b). 

Full sky maps provide the smallest record of the CMB anisotropy without loss of information. They 
permit a wide variety of statistics to be computed from the data - one of the most fundamental is the 
angular power spectrum of the CMB. Indeed, if the temperature fluctuations are Gaussian, with random 
phase, then the angular power spectrum provides a complete description of the statistical properties of the 
CMB. Komatsu et al. (2003) have analyzed the first-year WMAP sky maps to search for evidence of non- 
Gaussianity and find none, aside from a modest level of point source contamination which we account for in 
this paper. Thus, the measured power spectrum may be compared to predictions of cosmological models to 
develop constraints on model parameters. 

This paper presents the angular power spectrum obtained from the first-year WMAP sky maps. Com- 
panion papers present the maps and an overview of the basic results (Bennett et al. 2003b), and describe 
the foreground removal process that precedes the power spectrum analysis (Bennett et al. 2003c). Spergel 
ct al. (2003), Verde et al. (2003), Peiris et al. (2003), Kogut ct al. (2003), and Page ct al. (2003b) discuss the 
implications of the WMAP power spectrum for cosmological parameters and carry out a joint analysis of the 
WMAP spectrum together with other CMB data and data from large-scale structure probes. Hinshaw et al. 
(2003b), Jarosik et al. (2003b), Page et al. (2003a), Barnes et al. (2003), and Limon et al. (2003) discuss 
the data processing, the radiometer performance, the instrument beam characteristics and the spacecraft 
in-orbit performance, respectively. 

A sky map AT(n) defined over the full sky can be decomposed in spherical harmonics 



where n is a unit direction vector, and Yi m (n) are the spherical harmonic functions evaluated in the direction 
n. 




(1) 
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with 




(2) 
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If the CMB temperature fluctuation AT is Gaussian distributed, then each a; m is an independent 
Gaussian deviate with 

(ai m ) = 0, (3) 

and 

(ai m a*, m ,) = 6u> 5 mm > Q, (4) 

where Ci is the ensemble average power spectrum predicted by models, and S is the Kronecker symbol. The 
actual power spectrum realized in our sky is 

m— — I 

In the absence of noise, and with full sky coverage, the right hand side of equation (5) provides an unbiased 
estimate of the underlying theoretical power spectrum, which is limited only by cosmic variance. However, 
realistic CMB anisotropy measurements contain noise and other sources of error that cause the quadratic 
estimator in equation (5) to be biased. In addition, while WMAP measures the anisotropy over the full sky, 
the data near the Galactic plane are sufficiently contaminated by foreground emission that only a portion of 
the sky (~85%) can be used for CMB power spectrum estimation. Thus, the integral in equation (2) cannot 
be evaluated as such, and other methods must be found to estimate C/. 

In Appendix A, we review two methods that have appeared in the literature for estimating the angular 
power spectrum in the presence of instrument noise and sky cuts. The first (Hivon et al. 2002) is a quadratic 
estimator that evaluates equation (2) on the cut sky yielding a "pseudo power spectrum" C/. The ensemble 
average of this quantity is related to the true power spectrum, C; by means of a mode coupling matrix 
Gii' (Hauser & Peebles 1973). The second method (Oh et al. 1999) uses a maximum likelihood approach 
optimized for fast evaluation with WMAP-like data. In Appendix A we demonstrate that the two methods 
produce consistent results. 

A quadratic estimator offers the possibility of computing both an "auto-power" spectrum, proportional 
to J2 m \ a im\ 2 , and a "cross-power" spectrum, proportional to X^m a im a Zm' where the ai m coefficients are 
estimated from two independent CMB maps, i and j. This latter form has the advantage that, if the noise in 
the two maps is uncorrelated, the quadratic estimator is not biased by the noise. For all cosmological analyses, 
we use only the cross-power spectra between statistically independent channels. As a result, the angular 
spectra are, for all intents and purposes, independent of the noise properties of an individual radiometer. 
This is analogous to interferometric data which exhibits a high degree of immunity to systematic errors. The 
precise form of the estimator we use is given in Appendix A. 

The plan of this paper is as follows. In §2 we review the properties of the WMAP instrument and how 
they affect the derived power spectrum. In §3 we present results for the angular power spectra obtained from 
individual pairs of radiometers, the cross-power spectra, and examine numerous consistency tests. In §4, in 
preparation for generating a final combined power spectrum, we present the full covariance matrix of the 
cross-power spectra. In §5 we present the methodology used to produce the final combined power spectrum 
and its covariance matrix. In §6 we compare the WMAP first-year power spectrum to a compilation of 
previous CMB measurements and to a prediction based on a combination of previous CMB data and the 
2dFGRS data. We summarize our results in §7 and outline the power spectrum data products being made 
available through the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Appendix A 
reviews two methods used to estimate the angular power spectrum from CMB maps. Appendix B describes 
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how we account for point source contamination. Appendix C presents our approach to combining multi- 
channel data. Appendix D describes how the foreground mask correlates multipole moments in the Fisher 
matrix, and Appendix E collects some useful properties of the spherical harmonics. 

2. INSTRUMENTAL PROPERTIES 

The WMAP instrument is composed of 10 "differencing assemblies" (DAs) spanning 5 frequencies from 
23 to 94 GHz (Bennett et al. 2003a). The 2 lowest frequency bands (K and Ka) are primarily Galactic 
foreground monitors, while the 3 highest (Q, V, and W) are primarily cosmological bands. There are 8 
high frequency differencing assemblies: Ql, Q2, VI, V2, and Wl through W4. Each DA is formed from 
two differential radiometers which are sensitive to orthogonal linear polarization modes; the radiometers are 
designated 1 or 2 (e.g., Vll or W12) depending on which polarization mode is being sensed. 

The temperature measured on the sky is modified by the properties of the instrument. The most 
important properties that affect the angular power spectrum are finite resolution and instrument noise. Let 
Cf denote the auto or cross-power spectrum evaluated from two sky maps, i and i', where i is a DA index. 
Further, define the shorthand i = to denote a pair of indices, e.g., (Q1,V2). This spectrum will have 

the form 

C\ = w\C^ + Nl (6) 

where w\ = b\ b\ pf is the window function that describes the combined smoothing effects of the beam and 
the finite sky map pixel size. Here b\ is the beam transfer function for DA i, given by Page et al. (2003a) 
[note that they reserve the term "beam window function" for (b\) 2 ], and pi is the pixel transfer function 
supplied with the HEALPix package (Gorski et al. 1998). N} is the noise spectrum realized in this particular 
measurement. On average, the observed spectrum estimates the underlying power spectrum, Cj, 

(C\)=w\C l + {Ni)6 w , (7) 

where (Nf) is the average noise power spectrum for differencing assembly i, and the Kronecker symbol in- 
dicates that the noise is uncorrelated between differencing assemblies. To estimate the underlying power 
spectrum on the sky, C; , the effects of the noise bias and beam convolution must be removed. The determi- 
nation of transfer functions and noise properties are thus critical components of any CMB experiment. 

In §2.1 we summarize the results of Page et al. (2003a) on the WMAP window functions and their 
uncertainties. We propagate these uncertainties through to the final Fisher matrix for the angular power 
spectrum. In §2.2 we present a model of the WMAP noise properties appropriate to power spectrum 
evaluation. For cross power spectra (i ^ i' above), the noise bias term drops out of equation (7) if the noise 
between the two DAs is uncorrelated. These cross-power spectra provide a nearly optimal estimate of the 
true power spectrum, essentially independent of errors in the noise model, thus we use them exclusively in 
our final power spectrum estimate. The noise model is primarily used to propagate noise errors through 
the analysis, and to test a variety of different power spectrum estimates for consistency with the combined 
cross-power spectrum. 

2.1. Window Functions 

As discussed in Page et al. (2003a), the instrument beam response was mapped in flight using obser- 
vations of the planet Jupiter. The signal to noise ratio is such that the response, relative to the peak of 
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the beam, is measured to approximately —35 dB in W band, the band with the highest angular resolution. 
The beam widths, measured in flight, range from 0?82 at K band down to 0?20 in some of the W band 
channels (FWHM). Maps of the full two-dimensional beam response are presented in Page et al. (2003a), 
and are available with the WMAP first-year data release. The radial beam profiles obtained from these 
maps have been fit to a model consisting of a sum of Hermitc polynomials that accurately characterize the 
main Gaussian lobe and small deviations from it. The model profiles are then Legendre transformed to 
obtain the beam transfer functions b\ for each DA i. Full details of this procedure are presented in Page 
et al. (2003a), and the resulting transfer functions are also provided in the first-year data release. We have 
chosen to normalize the transfer function to 1 at I = 1 because WMAP calibrates its intensity response using 
the modulation of the CMB dipolc (I = 1). This effectively partitions calibration uncertainty from window 
function uncertainties. 

The beam processing described above provides a straightforward means of propagating the noise uncer- 
tainty directly from the time-ordered data through to the final transfer functions. The result is the covariance 
matrix S| u , for the normalized transfer function. Plots of the diagonal elements of Y,\ u , are presented in 
Page et al. (2003a). The fractional error in the transfer functions b\ are typically 1-2% in amplitude. In 
the end, these window function uncertainties dominate the small off-diagonal elements of the final covari- 
ance matrix for the combined power spectrum (see §4) . Additional observations of Jupiter will reduce these 
uncertainties. 

An additional source of error in our treatment of the beam response arises from non-circularity of the 
main beam. The effects of this non-circularity are mitigated by WMAP's scan strategy which results in most 
sky pixels being observed over a wide range of azimuth angles. The effective beam response on the sky is thus 
largely symmetrized. We estimate that the effects of imperfect symmetrization produce window function 
errors of < 1% relative to a perfectly symmetrized beam window function (Page et al. 2003a; Hinshaw et al. 
2003b). This error is well within the formal uncertainty given in In §5 we infer the optimal power 

spectrum by combining the 28 cross-power spectra measured by the 8 high frequency DAs Ql through W4. 
As part of this process we marginalize over the window function uncertainty, which automatically propagates 
these errors into the final covariance matrix for the combined power spectrum. Both the combined power 
spectrum and the corresponding Fisher matrix are part of the first-year data release. 

2.2. Instrument Noise Properties 

The noise bias term in equation (7) is the noise per ai m mode on the sky. If auto-power spectra are 
used in the final power spectrum estimate, the noise bias term must be known very accurately because it 
exponentially dominates the convolved power spectrum at high I. If only cross correlations are used, the 
noise bias is only required for propagating errors. Our final best spectrum is based only on cross correlations, 
and is independent of this term. However, as an independent check of our results, we evaluate the maximum 
likelihood spectrum based on a combined Q+V+W sky map. The noise bias must be estimated accurately 
for this application. 

In the limit that the time-ordered instrument noise is white, the noise bias is a constant, independent of I. 
If the time-ordered noise has a 1// component, the bias term will rise at low I. In this subsection we estimate 
the noise bias properties for each of the high frequency WMAP radiometers based on the time-ordered noise 
properties presented in Jarosik et al. (2003b). While the WMAP radiometer noise is nearly white by design 
(Jarosik et al. 2003a) with 1// knee frequencies of less than 10 mHz for 9 out of 10 differencing assemblies, 
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one of the radiometers (W41) has a 1// knee frequency of ~45 mHz. The latter is large enough that the 
deviations of (Nf) from a constant must be accounted for. 

The most reliable way to estimate the effects of 1// noise on the measured power spectra is by Monte 
Carlo simulation. Using the pipeline simulator discussed in Hinshaw et al. (2003b) we have generated a 
library of noise maps with flight-like properties. Specifically we have included flight-like 1// noise in the 
simulated time-ordered data, and have run each full-year realization through the map-making pipeline, 
including the baseline pre- whitening discussed in Hinshaw et al. (2003b). We evaluate the power spectra of 
these maps using the quadratic estimator described in Appendix A with 3 different pixel weighting schemes. 
(See Appendix A. 1.2 for definitions of the weights, and the I range in which each is used.) We define the 
effective noise as a function of I based on fits to these Monte Carlo noise spectra. For the analyses in this 
paper, we fit the spectra to a model of the form 



where the c l n are fit coefficients given in Table 1, with n max = 3 for I < 200, and n max = 1 for I > 200. 

Figure 1 shows the noise spectrum derived from the simulations for each of the 8 high frequency DAs, 
using uniform weighting over the entire I range. For comparison, we also plot an estimate of the CMB power 
spectrum from §5 in grey. Note that the W4 spectrum is the only one of this set to exhibit deviations from 
white noise in an I range where the signal-to-noise is relatively low, and we believe this simulation slightly 
over-estimates the 1/f noise in the flight W4 differencing assembly (Hinshaw et al. 2003a). 



Hinshaw et al. (2003a) present limits on systematic errors in the first-year sky maps. They consider the 
effects of absolute and relative calibration errors, artifacts induced by environmental disturbances (thermal 
and electrical), errors from the map-making process, pointing errors, and other miscellaneous effects. The 
combined errors due to relative calibration errors, environmental effects, and map-making errors are limited 
to < 15 /iK 2 (2cr) in the quadrupole moment C2 in any of the 8 high-frequency DAs. Tighter limits are 
placed on higher-order moments. We conservatively estimate the absolute calibration uncertainty in the 
first-year WMAP data to be 0.5%. 

Random pointing errors are accounted for in the beam mapping procedure; the beam transfer functions 
presented by Page et al. (2003a) incorporate random pointing errors automatically. A systematic pointing 
error of ~10" at the spin period is suspected in the quaternion solution that defines the spacecraft pointing. 
This is much smaller than the smallest beam width (~12' at W band), and we estimate that it would produce 
<1% error in the angular power spectrum at I = 1000, thus we do not attempt to correct for this effect. 
Barnes et al. (2003) place limits on spurious contributions due to stray light pickup through the far sidelobes 
of the instrument. They place limits of <10 ^K 2 on spurious contributions to Cj, at Q through W band, 
due to far sidelobe pickup. 

A detailed model of Galactic foreground emission based on the first-year WMAP data is presented by 
Bennett et al. (2003c) and is summarized in §3.1. We show that diffuse foreground emission is a modest 
source of contamination at large angular scales (> 2°). Systematic errors on these angular scales are negligible 
compared to the (modest) level of foreground emission. On smaller angular scale (< 2°), the 1-3% uncertainty 
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2.3. Systematic Errors 
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in the individual beam transfer functions is the largest source of uncertainty, while for multipolc moments 
greater than ^600, random white noise from the instrument is the largest source of uncertainty. 

3. THE DATA 

Figure 2 shows the cross-power spectra obtained from all 28 combinations of the 8 differencing assemblies 
Ql through W4 using the quadratic estimator described in Appendix A.l. These spectra have been evaluated 
with the Kp2 sky cut described in Bennett ct al. (2003c). The spectra are color coded by effective frequency 
V v i v %l ', where v % is the frequency of differencing assembly i. The low frequency (41 GHz) data are shown 
in red, the high frequency (94 GHz) data in blue, with intermediate frequencies following the colors of the 
rainbow. The top panel shows 1(1 + l)Ci/2n in /iK 2 , while the bottom panel plots the ratio of each channel 
to our final combined spectrum presented in §5. The top panel shows a very robust measurement of the first 
acoustic peak with a maximum near I ~ 220 and a shape that is consistent with the predictions of adiabatic 
fluctuation models. There is also a clear indication of the rise to a second peak at I ~ 540. See Page et al. 
(2003b) for an analysis and discussion of the peaks and troughs in the first-year WMAP power spectrum. 

The red data in the top panel show very clearly that the low frequency data are contaminated by diffuse 
Galactic emission at low / and by point sources at higher I. The higher frequency data show less apparent 
contamination, consistent with the foreground emission being dominated by radio emission, rather than 
thermal dust emission, as expected in this frequency range. 

3.1. Galactic and Extragalactic Foregrounds 

Bennett ct al. (2003c) present a detailed model of the Galactic foreground emission based on a Maximum 
Entropy analysis of all 5 WMAP frequency bands, in combination with external tracer templates. They 
demonstrate that the emission is well modeled by three distinct emission components. 1) Synchrotron 
emission from cosmic ray electrons, with a steeply falling spectrum in the WMAP frequency range: Ta(v) cx 
with < — 3, steepening with increasing frequency. 2) Free- free emission from the ionized interstellar 
medium that is well traced by Ha emission in regions where the dust extinction is low. 3) Thermal emission 
from interstellar dust grains with an emissivity index <~ 2.2. The model has a Galactic signal minimum 
between V and W band. 

In principal we could subtract the above model from each WMAP channel and recompute the power 
spectrum. However, since the model is based on WMAP data that have been smoothed to an angular 
resolution of 1?0, the resulting maps would have complicated noise properties. For the purposes of power 
spectrum analysis, we adopt a more straightforward approach based on fitting foreground tracer templates 
to the Q, V, and W band data. The details of this procedure, the resulting fit coefficients, and a comparison 
of the fits to the Maximum Entropy model are given in Bennett et al. (2003c). They estimate that the 
template model removes ~85% of the foreground emission in Q, V, and W bands and that the remaining 
emission constitutes less than ~2% of the CMB variance (up to I = 200) in Q band, and 1 1% of the CMB 
variance in V and W bands. 

The contribution from extragalactic radio sources has been analyzed in three separate ways. Bennett 
et al. (2003c) directly fit for sources in the sky maps. The result of this analysis is that we have identified 208 
sources in the WMAP data with sufficient signal to noise ratio to pass the detection criterion (we estimate 
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that ~5 of these are likely to be spurious). The derived source count law is consistent with the following 
model for the power spectrum of the unresolved sources 




with A = 0.015 /xK 2 sr (measured in thermodynamic temperature), (5 = —2.0, and vq = 45 GHz. Komatsu 
ct al. (2003) evaluate the bispectrum of the WMAP data and are able to fit a non-Gaussian source component 
based on a particular configuration of the bispectrum. They find the same source model, equation (9), fits 
the bispectrum data. For the remainder of this section, we adopt this model for correcting the cross- 
power spectra. At Q band (41 GHz) the correction to 1(1 + 1)Ci/2tt is 868 and 3468 /xK 2 at Z=500 and 
1000, respectively. At W band (94 GHz), the correction is only 31 and 126 /xK 2 at the same / values. 
For comparison, the CMB power in this I range is ~ 2000 /iK 2 . Later, when we derive a final combined 
spectrum from the multi-frequency data, we adopt equation (9) as a model with A as a free parameter. We 
simultaneously fit for a combined CMB spectrum and source amplitude and marginalize over the residual 
uncertainty in A. The best-fit source amplitude from this process is consistent with the other two methods. 

Figure 3 shows the cross-power spectra obtained from the same 28 combinations as in Figure 2, this 
time with the Galactic template model and source model subtracted. The bottom panel of the Figure shows 
the ratio of the 28 channels to the combined spectrum obtained in §5. The 28 cross-power spectra are 
consistent with each other at the 5 to 20% level over the I range 2 - 500. Similar scatter is seen in Monte 
Carlo simulations of an ensemble of 28 cross-power spectra with WMAP's beam and noise properties. The 
only significant deviation lies in the Q band data at low I which is ~40% below the higher frequency bands 
at I < 20. This is consistent with the accuracy estimated above for the Galactic template model, see also 
Figure 11 of Bennett et al. (2003c) for images of the maps after Galactic template subtraction. Since the 
WMAP data are not noise limited at low I, we use only V and W band data in the final combined spectrum 
for I < 100. 

The subtraction of the source model, equation (9), brings the Q band spectrum into good agreement 
with the other cross-power spectra up to I ~ 500. At higher I, the Q band data contributes very little to 
the final combined spectrum because the (normalized) Q band window function has dropped to less than 
5% (Page et al. 2003a). As discussed in §5, we marginalize over the source amplitude uncertainty, SA, when 
obtaining the final power spectrum estimate and associated covariance matrix. Thus the uncertainty is also 
accounted for in subsequent cosmological parameter fits (Verde et al. 2003; Spcrgel et al. 2003; Peiris et al. 
2003). 

Figure 4 shows a close-up of the 28 cross-power spectra in Figure 3 up to I — 100. The top panel shows 
the raw (unbinned) data which has correlations of <2% between neighboring points is this I range (see §4). 
These data are strikingly consistent with each other and support the conclusion that systematic errors at 
low I arc insignificant. To assess the level of scatter that does exist between the spectra, we have generated 
a Monte Carlo simulation in which we compute the rms scatter among the 28 spectra at each /, relative to 
the measured power. The bottom panel of Figure 4 shows the results of this simulation, averaged over 1000 
realizations, compared to the relative rms scatter in the data. The agreement is excellent, indicating that 
the uncertainty in the measured power spectrum in this I range is a few percent and is consistent with a 
combination of instrument noise and mode coupling due to the 15% sky cut. 

Another striking feature is the low amplitude of the observed quadrupole, and the sharp rise in power, 
almost linear in I, to I = 5. Bennett et al. (2003b) quote a value for the rms quadrupole amplitude, 
Q rms = \J (5/47r)C 2 = 8 ± 2 /xK, where the uncertainty is largely due to Galactic model uncertainty. This is 
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consistent with the amplitude measured by the COBE-DMR experiment, Q rms — W_ 4 /zK (Bennett et al. 
1996). The fast, nearly linear rise to I = 5 produces an angular correlation function with essentially no power 
on angular scales >60°, again in excellent agreement with the COBE-DMR correlation function (Bennett 
et al. 2003b; Hinshaw et al. 1996). In the context of a standard ACDM model, the probability of observing 
this little power on scales greater than 60° is ~ 2 x 10~' ! (Bennett et al. 2003b). 



4. THE FULL COVARIANCE MATRIX 

In §5 we derive our best estimate of the angular power spectrum by optimally combining the 28 cross- 
power spectra discussed above. The procedure for combining spectra requires the full covariance matrix of 
the individual cross-power spectra - in this section we outline the salient features of this matrix. There are six 
principal sources of variance for the measured spectra, C\: cosmic variance, instrument noise, mode coupling 
due to the foreground mask, point source subtraction errors, uncertainty in the beam window functions, and 
an overall calibration uncertainty. We ignore uncertainties in the diffuse foreground correction since they 
are everywhere sub-dominant to the cosmic variance uncertainty (see §3.1). 

We may write the covariance matrix as 

(£fuii)!i - ([C\ - (Q + A5>j][4 - (CV + A&)4 ]) , (10) 

where the angle brackets represent an ensemble average, C/ is the true underlying power spectrum, w\ is 
the window function of spectrum i, and AS 1 is the point source contribution to spectrum i. Here we have 
defined a point source spectral function, S 1 , as 



where z/j, vq, and (3 are as defined after equation (9). Note that (£f u ii)J;/ is symmetric in both {IV) and (ij). 

In the process of forming the combined spectrum we will estimate a best-fit point source amplitude, A, 
and subtract the corresponding source contribution from each spectrum i. We thus rewrite Ef u n as 

(S full )ii - (\C\ - (Q + SAS^wUlGi - {C V + SAS>)4 ]) , (12) 

where C\ = C\ — AS*w\ is the source-subtracted spectrum, and 8 A = A — A is the residual source amplitude, 
which we will marginalize over as a nuisance parameter. 

We may expand the covariance matrix as 

Efull ^cv ^ ^mask ^sre ^win- 

(13) 

We discuss each of these contributions in more detail below. 

Cosmic Variance, Noise, and Mode Coupling. — The first two terms, E cv + S mas k, incorporate the 
combined uncertainty due to cosmic variance, instrument noise, and mode coupling due to the foreground 
mask, 

(E cv + £ mask )!j, = ( [C\ - Q w\ ] [Cj, - C V 4 ] ) , (14) 

where w) is fixed at its measured value. We have split this contribution into two pieces to mimic the procedure 
we actually use to compute the covariance matrix. As outlined in §5.1, we start with S cv , then incorporate 
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the effects of point source error and window function error. We do not add the effects of mode coupling, 
Smask, until the very end of the computation. We consider this term in more detail in Appendix D. 

Point Source Subtraction Errors. — The third term, S src , is due to uncertainty in the point source 
amplitude determination 

(£ src )j), = S i w\al c 4,S\ (15) 

where cr 2 rc = ((A — A) 2 ) is the variance in the best-fit amplitude A, and we assume that the frequency 
dependence, S 1 , is perfectly known. In practice, we do not explicitly evaluate S src as given above, rather we 
employ a method based on marginalizing a Gaussian likelihood function, £(C/, A), over a nuisance parameter 
A. This process, which is discussed in Appendix B, effectively yields S cv + S src . 

Window Function and Calibration Errors. — The fourth term, S w ; n , is due to uncertainty in the beam 
window function, w\. This term arises from fluctuations in the window function which cause the measured 
spectrum, C\ to differ from our estimate of the convolved spectrum, CjioJ, where w\ is the estimated window 
function. This contribution has the form 

(S win )« = C x (At«}.Ai4) C v . (16) 

Recall from §2 that w\ — b\b\ pf where b\ is the beam transfer function for DA i and pf is the pixel window 
function. Define u\ = Abj/bj to be the fractional error in b\, then to first order in u\ we have 

(Aw\ ■ A4 ) = 4 + + 4 4 + 44 ) 4> • ( 17 ) 

For WMAP the uncertainty in b\ is uncorrelated between DAs, thus the above expression reduces to 

■ a4) = w \ [k,u'(^ + + KwQi'i + s i'j')} 4 > (is) 

where we have defined 

Ei iIJ , = («?4) = (6}6j,)- 1 (19) 

and where u , is the covariance matrix of the beam transfer function for DA i given by Page et al. (2003a). 
When generating the combined spectrum in the next section, we add S w ; n to the above contributions, giving 

S cv + S src -(- Xj w j n . 

The WMAP absolute calibration uncertainty is 0.5%. We do not explicitly incorporate this contribution 
in the covariance matrix. Instead, we propagate a 0.5% uncertainty into the normalization of the final power 
spectrum amplitude (Spergel et al. 2003). 



5. THE COMBINED POWER SPECTRUM 

In §3 we demonstrated that the three high frequency bands of WMAP data produced consistent estimates 
of the angular power spectrum, after a modest correction for diffuse Galactic emission and extragalactic point 
sources. It is therefore justifiable to combine these data into a single "optimal" estimate of the angular power 
spectrum of the CMB. In this section, we provide an overview of two methods we use to generate a single 
combined spectrum. The first is a multi-step process that simultaneously fits the 28 cross-power spectra 
presented above to a single CMB power spectrum and a point source model, equation (9), while correctly 
propagating beam and residual point source uncertainties through to a final Fisher matrix. This spectrum 
constitutes our best estimate of the CMB power spectrum from the first-year WMAP data. The second 
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spectrum, which serves as a cross check of the first, is based on forming a single co-added sky map from the 
Ql through W4 maps, and using the quadratic estimator with noise bias subtraction. We compare the two 
spectra in §5.3. 



5.1. Method I Optimal Combination of Cross-Power Spectra 

Since this method is relatively complicated, we outline the basic procedure here and relegate the details 
to Appendices, as indicated. We present the result in §5.3. The steps are as follows. 

1. Subtract best-fit Galactic foreground templates from each of the maps Ql through W4, using the 
coefficients given in Table 3 of Bennett et al. (2003c). 

2. Evaluate the 28 cross-power spectra from the maps Ql through W4, where each spectrum has been 
evaluated using the quadratic estimator of Appendix A.l with the weighting scheme defined in Ap- 
pendix A. 1.2. 

3. Collect the noise bias estimate, ( A ; l ) , for each DA from §2.2. These estimates are used in the calculation 
of the covariance matrix for the combined spectrum, and in setting the relative weight of each cross- 
power spectrum in the final combined spectrum. 

4. Apply the procedure presented in Appendix B.l to obtain an estimate for the point source amplitude. 
The value obtained is A — 0.0155 ± 0.0017, roughly independent of l max in equation (B6). This value 
is in good agreement with an estimate based on the bispectrum (Komatsu et al. 2003), and on an 
extrapolation of point source counts (Bennett et al. 2003b). Subtract the point source contribution 
from the cross-power spectra: C\ = C\ — AS l w\. 

5. Compute an approximate form of the full covariance matrix discussed in §4, £f u ii. The procedure we 
use produces a covariance matrix which includes cosmic variance, instrument noise, source subtraction 
uncertainties, and window function uncertainties. At this stage of the process, it does not yet include 
the effects of mode coupling. More details are given in §4 and Appendix B. 

6. Invert the approximate covariance matrix for use in computing the optimal spectrum. This is the most 
computationally intensive step in the process. 

7. Compute the final combined spectrum from the 28 C\ as per the procedure given in Appendix C. In 
particular, assume a fiducial cosmological model (as specified in the Appendix), and use equation (C6), 
with ;; , — w . This produces a final spectrum which is very nearly optimal. 

Note: for I < 100 we use a surrogate procedure for computing the combined spectrum. In order to 
minimize Galactic foreground contamination, we use only V and W band data. Moreover, because 
the statistics of the C\ are mildly non-Gaussian, and because point source contamination and window 
function uncertainties are small, the "optimal" machinery developed above is unnecessary Rather we 
simply form a weighted average spectrum from the V and W band C\. 

8. Compute the approximate inverse-covariance matrix, Q, for the combined spectrum using equation (C7). 
This matrix is approximate in two ways. 1) It docs not yet incorporate the effects of mode coupling 
- this is added below. 2) It has been evaluated for a fiducial cosmological model, Q(Cf d ) 1 while in a 
likelihood application, we need to evaluate <5(C ; th ) for an arbitrary model, Cf h - we add this next. 
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9. Introduce the dependence on cosmological model into Q as follows. Invert Q to obtain the approximate 
covariance matrix of the combined spectrum, S. The off-diagonal terms of £ are small and weakly 
dependent on cosmological model, so we expand S as 

t w =D l 5 w +e w , (20) 

where ew encodes the mode coupling due to window function and source subtraction uncertainties. 
This relation defines ew , which we take to be zero on the diagonal. Di is dominated by cosmic variance 
and noise; in order to separate the two contributions, we introduce an ansatz of the form 

Dl ^zrrij^ y {cfd+Nf)2 ' (21) 

where Cf d is the fiducial model used to generate the combined spectrum (see Appendix C), and Nf s 
is the effective noise in the combined spectrum, which is defined by this equation. The dependence 
on cosmological model is introduced in the covariance matrix by re-computing Di with Cf d — > Cf h , 
leaving Nf s fixed. 

10. Estimate the coupling induced by the foreground mask as described in Appendix D. The effect of the 
mask on the off-diagonal elements of the Fisher (inverse-covariance) matrix can be written as 

rwiDW,)- 1 '*, (22) 

where D[ = Dif~^. This expression parameterizes, and defines, the off-diagonal mode coupling in 
the form of a correlation matrix, rw ■ Note that Fw, defined in Appendix D.2, is related to rw by 
Fw = S w +r w . 

11. The final Fisher (curvature) matrix, Qw, is obtained using 

Qw = {D[)- l 5 lv e w (£>;£>£' ) 1/2 + rw (A'^)" 1/2 , (23) 

We further calibrate D[ with Monte Carlo simulations. This calibration process, and a description of 
how the curvature matrix is used in a maximum likelihood determination of cosmological parameters, is given 
in §2 of Verde et al. (2003). As part of the first-year data release, we provide a Fortran 90 subroutine to 
evaluate the likelihood of a given cosmological model, Cf h , given the WMAP data (supplied in the routine). 
The code also optionally returns the Fisher (inverse-covariance) matrix for the combined spectrum. 



5.2. Method II Combined Sky Map 

Our second approach is to form a single co- added map from the Ql through W4 maps, 

V^lO rpl I 2 

t = t:r 3 ;; r > (24) 

where T[ is the sky map for DA i with the best-fit Galactic template model subtracted, and (Tq i is the noise 
per observation for DA i, given by Bennett et al. (2003b). We evaluate the power spectrum of this map 
on the Kp2 cut sky using the quadratic estimator in Appendix A.l. An effective noise model is obtained 
by using the same approach as described in §2.2: we generate co- added noise maps from the library of 
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end-to-end simulations, evaluate their average spectra, then fit a noise model. The noise bias model is then 
subtracted from the power spectrum of the combined temperature map. We have performed this analysis 
with 3 distinct pixel weighting schemes (see Appendix A. 1.2) and 3 corresponding noise models. The results 
are shown in Figure 5 where it is seen that the three cases are virtually indistinguishable. 

In effect, this analysis uses both the auto- and cross-power spectra. We view this as a useful check of the 
more sophisticated procedure described in §5.1, but we do not rely on it for a final result. Uncertainties in 
the noise model only effect the fourth moment of the cross-power spectra, but they effect the second moment 
of the auto-power spectra and potentially bias the final result. The ^6% sensitivity advantage gained by 
including auto-power spectra was not deemed worth the effort required to guarantee that the final result was 
not biased. 



5.3. Comparison of Results 

Figure 6 compares the power spectra obtained from methods I and II above. The combined cross-power 
spectrum from §5.1 is shown in black, the auto-power spectrum obtained in §5.2 from the co-added map is 
shown in grey. The two methods agree extremely well with the only notable deviation being at the highest I 
range probed by the first-year data. This is the regime where the auto-power spectrum will be most sensitive 
to the noise bias subtraction. As can be seen in the error estimates shown in Figure 8, the deviation is less 
than lcr. 

A separate test of robustness is to compute the angular power spectrum in separate regions of the sky to 
see if the spectrum changes. We have computed the power spectrum in two subsets of the sky - the ecliptic 
poles, and the ecliptic plane, using the quadratic estimator with the combined sky map. The results are 
shown in Figure 7 where the pole data is shown in grey and the plane data in black. The two spectra are 
very consistent overall, but some of the features that appear in the combined spectrum, such as the "peak" 
at I ~ 40 and the "bite" at I ~ 210, are not robust to this test, thus we consider these features to be of 
marginal significance. There is also no evidence that beam ellipticity, which would be more manifest in the 
plane than in the poles, systematically biases the spectrum. This is consistent with estimates of the effect 
given by Page et al. (2003a). 

6. DISCUSSION 

Our best estimate of the angular power spectrum of the CMB is shown in Figure 8. Also shown is the 
best-fit ACDM model from Spergel et al. (2003) which is based on a fit to the this spectrum plus a compilation 
of additional CMB and large-scale structure data. The WMAP data points are plotted with measurement 
errors based on the diagonal elements of the Fisher matrix presented in Appendix D. The cosmic variance 
errors, which include the effects of the sky cut, are plotted as a lcr band around the best-fit model. As 
discussed in Spergel et al. (2003), the model is an excellent fit to the data. The combined spectrum provides 
a definitive measurement of the CMB power spectrum with uncertainties limited by cosmic variance up to 
I <~ 350. The spectrum clearly exhibits a first acoustic peak at I = 220.1 ± 0.8 and a second acoustic peak 
at I = 546 ± 10. Page et al. (2003b) present an analysis and interpretation of the peaks and troughs in the 
first-year WMAP power spectrum. 

Figure 9 compares the first-year WMAP spectrum to a compilation of recent balloon and ground-based 
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measurements. In order to make this Figure meaningful, we plot the best-fit model spectrum to represent 
the WMAP results. The data points are plotted with errors that include both measurement uncertainty 
and cosmic variance, so no error band is included with the model curve. (Since individual groups report 
band power measurement with different band-widths, it is not possible to represent a single cosmic variance 
band that applies to all data sets.) The model spectrum fit to WMAP agrees very well with the ensemble 
of previous observations. 

Wang et al. (2002a) have recently distilled a CMB power spectrum from an optimal combination of 
the extent pre- WMAP data. In Figure 10 we plot their derived band power points alongside the WMAP 
data. To make this comparison meaningful, we plot the WMAP data with cosmic variance plus measurement 
errors and omit the error band from the model spectrum. The distilled spectrum is notably lower than the 
WMAP data in the vicinity of the first acoustic peak. In a previous version of this work (Wang et al. 2002b) 
the authors noted that the first peak of their combined spectrum was lower than a significant fraction 
of their input data. They attribute this to their formalism allowing for a renormalization of individual 
experiments within their respective calibration uncertainties. Figure 1 in Bennett et al. (2003a) presents a 
similarly distilled spectrum from the data extent in late 2001 and found a first peak amplitude that was 
more intuitively consistent with the bulk of the input data, and which is now seen to be consistent with the 
WMAP power spectrum. 

Figure 11 shows the WMAP combined power spectrum compared to the locus of predicted spectra, in 
red, based on a joint analysis of pre- WMAP CMB data and 2dFGRS large-scale structure data (Percival 
et al. 2002). As in Figure 8, the WMAP data are plotted with measurement uncertainties, and the best- fit 
ACDM model (Spergel et al. 2003) is plotted with a la cosmic variance error band. Percival et al. (2002) 
predict the location of the first peak should occur at I — 221.8 ± 2.4, which is quite consistent with the 
value reported by Page et al. (2003b) of / = 220.1 ± 0.8. The height of the first peak was predicted to 
be in the range 4920 ± 170 ^K 2 , while Page et al. (2003b) report a measured height of 5580 ± 75 /iK 2 , 
about 13% higher. Unlike the position, the amplitude of the first peak has a complicated dependence on 
cosmological parameters. Percival et al. (2002) report best-fit parameters for a ACDM model that are 
mostly consistent with those reported by Spergel et al. (2003) for the same class of models. The only mildly 
disparate comparison lies in the combination of normalization, <r s , and optical depth, r. Percival et al. (2002) 
report the product <T$e~ T = 0.72 ± 0.03 ± 0.02, where the first error is a "theory" error and the second is 
measurement error. While Spergel et al. (2003) does not report a maximum likelihood range for this explicit 
parameter combination, the product of their maximum likelihood values for erg and r yields crge -1 " — 0.74, 
which is consistent with Percival et al. (2002), but would make the first peak a few percent higher. Small 
differences in VLbh 2 , fl m h 2 , and n s , may also contribute to the difference. 



7. CONCLUSIONS 

We present measurements of the angular power spectrum of the cosmic microwave background from the 
first-year WMAP data. The eight high-frequency sky maps from DAs Ql through W4 were used to estimate 
28 cross-power spectra, which are largely independent of the noise properties of the experiment. These data 
were tested for consistency in §3, then used in §5 as input to a final combined spectrum, discussed in §6. 
The procedure for estimating the uncertainties in the final combined spectrum were discussed in §4 and in 
numerous Appendices. 

The combined spectrum provides a definitive measurement of the CMB power spectrum, with uncer- 
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taintics limited by cosmic variance up to I <~ 350, and a signal to noise per mode > 1 up to I <~ 650. The 
spectrum clearly exhibits a first acoustic peak at I = 220.1 ± 0.8 and a second acoustic peak at I = 546 ± 10. 
Page et al. (2003b) present an analysis and interpretation of the peaks and troughs in the first-year WMAP 
power spectrum. Spergel et al. (2003), Verde et al. (2003), and Peiris et al. (2003) analyze the combined spec- 
trum in the context of cosmological models. They conclude that the data provide strong support adiabatic 
initial conditions, and they give precise measurements of a number of cosmological parameters. Kogut et al. 
(2003) analyze the correlation between WMAP J s temperature and polarization signals, the Cf E spectrum, 
and present evidence for a relatively high optical depth, and an early period of cosmic reionization. Among 
other things, this result implies that the temperature power spectrum is suppressed by ~30% on degree 
angular scales, due to secondary scattering. 

A variety of first-year WMAP data products are being made available by NASA's new Legacy Archive 
for Microwave Background Data Analysis (LAMBDA). In addition to the sky maps and calibrated time- 
ordered data, we are providing the 28 cross power spectra used in this paper (with diffuse foregrounds 
subtracted), the combined spectrum from §5.1, and a Fortran 90 subroutine to compute the likelihood of a 
given cosmological model, (the code will also optionally return the Fisher (inverse-covariance) matrix for the 
combined spectrum.) The LAMBDA URL is http://lambda.gsfc.nasa.gov/. 

The WMAP mission is made possible by the support of the Office of Space Sciences at NASA Head- 
quarters and by the hard and capable work of scores of scientists, engineers, technicians, machinists, data 
analysts, budget analysts, managers, administrative staff, and reviewers. We thank Mike Nolta for helpful 
comments on an earlier draft of this manuscript. LV is supported by NASA through Chandra Fellowship 
PF2-30022 issued by the Chandra X-ray Observatory center, which is operated by the Smithsonian Astro- 
physical Observatory for an on behalf of NASA under contract NAS8-39073. We acknowledge use of the 
HEALPix package. 
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A. POWER SPECTRUM ESTIMATION METHODS 

For the analysis of WMAP's first-year data, we have chosen two distinct methods for inferring the power 
spectrum. The first is a fast and accurate quadratic method for estimating the power spectrum of a partial 
sky map (Hivon ct al. 2002). We summarize the basic approach here, highlighting the aspects of the method 
that are especially pertinent to WMAP, and refer the reader to Hivon et al. (2002) for details. The second is 
a maximum likelihood method that provides an independent estimate of the spectrum measured by WMAP 
(Oh et al. 1999). 

We have applied both of these methods to the WMAP data. The results are shown in Figure 12, 
which shows spectra estimated from the V band map, up to I — 200, for the two methods. The maximum 
likelihood estimate has slightly lower uncertainties at low I, because the method optimally weights the data 
with a pixel-pixel covariance C = S + N w S where S is the covariance of the CMB signal and N is the 
covariance of the noise (see Appendix A. 2). Our quadratic estimator uses uniform pixel weights at low I 
(see Appendix A. 1.2) though it is clear from the Figure that the difference is not significant. At high I, 
where the WMAP data are noise dominated, the two estimators give essentially identical results because 
they effectively weight the data in the same way. 

To obtain our "best" estimate of the WMAP power spectrum, we adopt the quadratic estimator because 
it can be readily applied to pairs of WMAP radiometers in a way that is nearly independent of the properties 
of the instrument noise. In §5 we discuss our methodology for combining spectra from pairs of radiometers 
and present the final combined spectrum. 



A.l. Quadratic Estimation 

Hivon et al. (2002) start with the full-sky estimator in equation (2), add a position dependent weight, 
W(n), and set W to zero in the regions where the sky is contaminated. In other regions, W is chosen to 
optimize the sensitivity of the estimator. A temperature map AT(n) on which a weight IT(n) is applied can 
be decomposed in spherical harmonic coefficients as 

ai m = ydO n AT(n)IT(n)r ; :„(n) (Al) 

« O^AT(„)lf( P )C( P ), (A2) 
p 

where the integral over the sky is approximated by a discrete sum over map pixels, each of which subtends 
solid angle f2 p . Hivon et al. (2002) then define the "pseudo power spectrum" C\ as 

m— — I 

The pseudo power spectrum Ci, given by the weighted spherical harmonic transform of a map, is clearly 
different from the full sky angular spectrum, C ; sky , but the ensemble averages of the two spectra can be 
related by 

(C l ) = Y,Gu,(C^) (A4) 
i' 

where Gw describes the mode coupling resulting from the weight function W(n) (Hauser & Peebles 1973). 
Hivon et al. (2002) give the following expression for the coupling matrix, which depends only on the geometry 
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of the weight function W(n) 

^^E^ + iW.fo' o o) 2 ' (A5) 

where the final term in parentheses is the Wigner 3-j symbol, and Wi is the angular power spectrum of the 
weight function 

m 

where 

wi m 

Upon inverting the coupling matrix Gw and making the identification (C ; sky ) — Ci, we obtain the 
following estimator of the power spectrum 

Q obs = E G «<CV. (A8) 
v 

The computation of equation (A2) for each (I, m) up to I = Z max would scale as A p i x Z^ ax if performed on an 
arbitrary pixelization of the sphere, where N p ^ is the number of sky map pixels. However, for a pixelization 
scheme with iso-latitude pixel centers, fast FFT methods may be employed to speed up the evaluation, 
so it scales like N^l^^ (Muciaccia et al. 1997). The WMAP sky maps have been produced using the 
HEALPix layout (Gorski et al. 1998) which supports such fast spherical harmonic transforms. In particular, 
the HEALPix routine map2alm evaluates equation (A2). 

A. 1.1. Auto- and Cross-Power Spectra from the WMAP Data 

For a multi-channel experiment like WMAP it is quite powerful to evaluate the power spectra from 
different maps and compare results. In particular, the quadratic estimator described above may be used on 
1 or 2 maps at a time by generalizing the expression for the pseudo power spectrum equation (A3) as 

^ = 27TT S (A9) 

m— — I 

where a\ m refers to the transform of map i and aj^ refers to map j, which needn't be the same as map i. As 
discussed in §2, if i ^ j and the noise in the two maps is uncorrelated, the estimator equation (A9) provides 
an unbiased estimate of the underlying power spectrum. 

We have tested the auto- and cross-power estimator extensively with Monte Carlo simulations of the 
first-year WMAP data. Selected results from this testing are shown in Figure 13. The auto- and cross-power 
spectra that obtain from the flight data are presented in detail in §3. The cross-power spectra form the basis 
for our final combined spectrum, presented in §5. 

A. 1.2. Choice of Weighting 

We seek a weighting scheme that mimics the maximum likelihood estimation (Appendix A. 2), which 
effectively weights the data by the full inverse covariance matrix C _1 = (S + N) -1 . For the combined 



\Wlm\ 2 , 



21 + 1^ 

m 

J dn n W{n)Y; m {n). 



(A6) 
(A7) 
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spectrum presented in §5, we use three distinct weighting functions in three separate I ranges. 

1. For I < 200 we give equal weight to all un-cut pixels, 

W(p) = M(p) (A10) 

where M(p) is the Kp2 sky mask, defined by Bennett et al. (2003c). It takes values of within the 
mask and 1 otherwise. 

2. For I > 450 we use inverse- noise weighting, 

W(p) = M(p)N ohs (p) (All) 
where N b s (p) is the number of observations of pixel p. 

3. For 200 < I < 450 we use a transitional weighting, 

w<r) " U(nJ?UKM (A12) 

where (-/V b s ) is the mean of A b s evaluated over the cut sky. 
Verde et al. (2003) discuss the choice of weighting in more detail. 

A. 2. Maximum Likelihood Estimation 

This Appendix provides a summary of the maximum likelihood approach to power spectrum estimation 
originally presented by Oh et al. (1999). If the temperature fluctuations are Gaussian, and the a priori 
probability of a given set of cosmological parameters is uniform, then the power spectrum may be estimated 
by maximizing the multi-variate Gaussian likelihood function 

nn\ \ exp(-4m T CT 1 m) 

r«7,|m) = 7 == (A13) 

where m is a data vector (see below) and C is the covariance matrix of the data, which has contributions from 
both the signal and the instrument noise, C — S + N. We can work in whatever basis is most convenient. 
In the pixel basis the data are the sky map pixel temperatures, and in the spherical harmonic basis the data 
are the ai m coefficients of the map. In the former basis the noise covariance is nearly diagonal, while in the 
latter, the signal covariance is. 

Pixel basis: Spherical harmonic basis: 

m — > Ti m — > ai m 

s - E,T^P,(cos%) S - diag(C 2 ,C 2 ,...C3,C 3 ,...) 

N -> af5 l3 N -> N {lm){lm y (sec below). 

For WMAP, the length of the data vector, iV data , is 2,672,361, the number of 7' sky map pixels (HEALPix 
Aside = 512) that survive the Galaxy cut, so it is necessary to find methods for evaluating C that do not 
require a full inversion of the covariance matrix C, which requires O (A| ata ) operations. We use an iterative 
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method for evaluating the likelihood that exploits the ability to find an approximate inverse C . The most 
important features in the data that make this possible are that WMAP observes the full sky and the Galaxy 
cut is predominantly azimuthally symmetric in Galactic coordinates (Bennett et al. 2003b). Of secondary 
importance for this pre-conditioner is that WMAP's noise per pixel does not vary strongly across the sky 
(Bennett et al. 2003a). We discuss the pre-conditioner in more detail below. 

Defining / = — 21n£ and P ; = J^f , we maximize the likelihood by solving 

^j- =0 = m T C- 1 P*Cr 1 ni + tr(C- 1 PO (A14) 
aCi 

using a Newton-Raphson root finding method that generates an iterative estimate of the angular power 
spectrum 

at each step. Here Fw is the Fisher matrix 

F w - - (d 2 InC/dddQ,) = ^(CT 1 P, CT 1 P,,). (A16) 
To implement the solution in equation (A15) we need a fast way to evaluate the following three components 

1. m T C X P ; C 

2. ixiC- 1 ?!) 

3. t^C^PiC- 1 ^). 

We use the spherical harmonic basis in which the data vector consists of the a; m coefficients of the map 
obtained by least squares fitting on the cut sky. The signal covariance is diagonal in this basis, while the 
noise matrix is obtained from the normal equations for the a/ m fit 

2 N (lm)(lmy a W = Wm) ( A17 ) 
(lm)' 

where 

1 _ \ - Y^jn^Y^yjn,) 

iV (im)(im)' = 2^ ^2 ( A18 ) 

i 1 

\- TiY {lm) {hi) 

i 1 

The sums are over all sky map pixels that survive the Galaxy cut, and we have assumed that the noise is 
uncorrelated from pixel to pixel. 



A. 2.1. Evaluation of C 1 m 

The term C _1 m appears repeatedly in the evaluation of equation (A15). We compute this by solving 
C z = (S + N) z = m for z. A more numerically tractable system is obtained by multiplying both sides by 
SjN" 1 , so 

(I + S'N- 1 S')S'z = S'N- 1 m = S*y (A20) 
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where y is the spherical harmonic transform of the map, defined in equation (A19). Note that y can be 
quickly computed in any pixelization scheme that has iso-latitude pixel centers with fixed longitude spacing. 
We then solve equation (A20) using an iterative conjugate gradient method with a pre-conditioner for the 
matrix A = (I + S2 N _1 S'). We find the following block-diagonal form of A to be a good starting point 

- _ / I + S^N^S* \ 

A ^ diagP + SiN^S*) ) (A21) 

where N _1 is a block-diagonal approximation of the noise matrix discussed below. The lower-right block 
of A occupies the high I portion of the matrix where the signal to noise ratio S^ N -1 S^ is low, so a 
diagonal approximation is adequate. The upper-left block occupies the low I portion of the matrix where 
the signal dominates the noise, so we need a better estimate of N^ 1 . In practice we find this split works 
well at I = 512 for the WMAP noise levels. As for the approximate form of N _1 , defined in equation (A19), 
note that the dominant off-diagonal contributions arise from the azimuthally symmetric Galaxy cut, which 
couples different I modes, but not m modes. Thus N -1 is approximately block diagonal, with perturbations 
induced by the non- uniform sky coverage of WMAP. We therefore use a block diagonal form of N _1 as the 
pre-conditioner, 

N(lm)(lm)' = N (lL)(lm)'^ mm ' '■ (A22) 

Using the pre-conditioner equation (A21) we find that the conjugate gradient solution of equation (A20) 
converges in approximately six iterations and requires only cpu-minutes of processing on an SGI Origin 
2000. 



B. POINT SOURCE SUBTRACTION 

In this Appendix we describe the procedure we use to estimate and subtract the point source contri- 
bution directly from the multi-frequency cross-power spectra. We then show how we incorporate the source 
model uncertainty into the covariance matrix of the source-subtracted spectra by marginalizing a Gaussian 
likelihood function over the source model amplitude parameter. 

This marks the starting point of the multi-frequency analysis which will lead to the combined power 
spectrum, discussed in §5.1. In order to generate the combined spectrum, we need the full covariance matrix 
of the cross-power spectra (see §4) . Our approach to generating the full covariance is to start with the ideal, 
full-sky form, which only includes cosmic variance and instrument noise, then we incorporate additional 
effects step by step, as outlined in §5.1 and in these Appendices. For an ideal experiment with full sky 
coverage, no point source contamination, and no beam uncertainty the covariance matrix is 

v"'jj' — ! 



(Bl) 



(2/ + 1) 

where Sij denotes the Kronecker symbol, wf = b\ b\ pf is the window function, and n l n l 6w = N 1 . 



B.l. Estimating the Point Source Amplitude 

We start by assuming a Gaussian likelihood for the sky model, given the WMAP data 

-2Inr(A,C,|C?)=2 [Cl-id+AS^w]] (XT 1 ^ [Cj - (Q + AS*)w\ }, (B2) 

ij I 
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where C\ is cross-power spectrum i, w\ is the window function for spectrum i, C; is the true CMB power 
spectrum, and AS" is the source model defined in equations (9) and (11). Here we assume the diagonal form 
of the covariance matrix in equation (Bl). 

To determine the best-fit source amplitude, A, we marginalize this likelihood over the CMB spectrum, 
C/. First we expand equation (B2) as 

-21n£ = 22(C/-^)(S- 1 )} 1 (qf-^) 

l ij 

- ^2C i ^(Q i -A5X)(S- 1 )j j W i 

l ij 

+ E C 'E ^(S" 1 )?^. ( B3 ) 
which is a quadratic form in C;, £ cx Jj ; exp[— ^(aCf — IbCi + c)], with 

ij 

b ee ^(Cj-A^^KE- 1 )^, 

ij 

c ee ^(Cj-A^^KS- 1 )) 1 ^-^). (B4) 



We wish to evaluate the marginalized likelihood, Cci (A) = j dCiC(A,Ci). This is readily evaluated using 
the substitution Ci — > Ci — b/a, giving £c t oc IL ex P[ — \{°~ b 2 /a)], where we drop a multiplicative term 
proportional to a, which is independent of A. The marginalized likelihood function is thus 

-21n£ C( = ^^Cf-AS^XE-^ql-AS^) 



E 



Ei'j' w \' (s- 1 );'- 1 '^' 

Setting d£ci (A) jdA = gives the most likely point sources amplitude 

A _ ZiiiCjp-^Hiwl 



^(Q-A^X)^- 1 )^ 



where 



Hi = S 



The standard error on the best fit value of A is 



(B5) 



(B6) 



(B7) 



(B8) 



For WMAP, the off-diagonal terms in the covariance matrix are small. Here, neglecting them changes the 
inferred point source amplitude by less than 1%. 
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B.2. Marginalizing over Point Source Amplitude 

The source subtraction procedure discussed above is uncertain. In this section we incorporate this un- 
certainty into the full covariance matrix for the cross-power spectra. We again assume a Gaussian likelihood 
function of the form 

-2\nC=Y [C\-{^ + 8AS\)w\] (XT 1 )^, [& v - (C?, h + 5AS\, )w{, } (B9) 
y w 

where C\ is the source-subtracted cross-power spectrum for DA pair i, obtained above, w\ is the window 
function for spectrum i, Cf h is now a fixed CMB model spectrum, and 8 A is the residual source amplitude. 

We can marginalize the likelihood function over the residual point source amplitude as follows. Expand 
equation (B9) as 

-21n£ = Y (Cl-C^iZ-^iCi-Cpvi,) 

ij W 

ij W 

+ (5AfJ2 S'wKE-^S^,, (BIO) 

ij IV 

which is a quadratic form in 6 A, C cx exp[— ^(a(SA) 2 — 2b(SA) + c)], with 

ij W 

b ee J2 (ci-crwi)(^yi,si4, 

ij W 

c ee (Cl-CrwD^-^iCl-C^,). (Bll) 

ij IV 

We wish to evaluate the marginalized likelihood, £a = J Cd{5A). This is readily evaluated using the 
substitution SA — > SA—b/a, giving La °c exp[— \(c— b 2 /a)), where we drop a multiplicative term proportional 
to a, which is independent of C\. The marginalized likelihood function is thus 

-21n£.4 = Yl (Cl-^wl^-^iCf,-^,) (B12) 

ij W 

- \Y {C\-C^w\){H-%S^ v Y ^i" ~ C r wi i" ){Z- 1 )%,,S y w 3 l ',, l . (B13) 

ij IV i'j' l"V" 

This expression may be recast in the form 

-2\nC A = Y (Cl-Ct h w\) (F^)% (Ci-C?4), (B14) 
ij w 

where (F src )% is 

(F^)% ee (^)% - \b% , (B15) 

with 

b%= y (s- 1 )!!"^'^;,^- 1 )^,,^'^:. (Bi6) 

i'j' 



:/:/ /////// 
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The superscript "src" indicates that the Fisher matrix so obtained includes point source subtraction un- 
certainty, in addition to whatever effects have been included in (E -1 )^, already, in this case only cosmic 
variance and noise. Note that equation (B15) neglects a term proportional to lndet a which has a weak 
dependence on cosmological parameters. In the actual calculation, as in the previous section, we assume the 
diagonal form of the covariance matrix in equation (Bl). 



C. OPTIMAL WEIGHTING OF MULTI-CHANNEL SPECTRA 

We use the 8 high-frequency differencing assemblies Ql through W4 to generate the final combined 
spectrum. In this Appendix we show how we combine these spectra into a single estimate of the angular 
power spectrum. 

The ultimate goal of the WMAP analysis is to produce a likelihood function for a set of cosmological 
parameters, {a}, given the data, C\. Specifically 

V{a\C\) = C{C\\Ct\a))V{$), (CI) 

where V{a\C\) is the probability of {a} given the data, £(Cj |Cj th (a)) is the likelihood of the data given the 
model, Cf h (a), and V(a) is the prior probability of the parameter set (Verde et al. 2003). To this end, we 
seek a combined spectrum, Ci, that estimates the power spectrum in our sky, Cf ky , with the property that 
T(d\6i) =V(d\C\), and hence C(Ci\Cf h (d)) = L{C\\C\ h {d)). 

To estimate the combined spectrum, we approximate the likelihood function for the cross-power spectra 
as Gaussian 

-21n£(C?|Cf) = £(C? - C?w\) (S^)« {C{, - Cf v? v ), (C2) 

ij W 

where C\ is the spectrum with the best-fit source model subtracted, defined after equation (12), w\ is the 
window function of spectrum i, defined after equation (6), and £f u ii is the covariance matrix of the 28 cross- 
power spectra. Note that the treatment in the remainder of this section does not depend on any specific 
property of the covariance, so we use generic notation for readability. However, when we generate the 
WMAP first-year combined spectrum, the actual form of the covariance used at this step is (Sf u ii)JJ/ , where 
£ indicates the approximate covariance, obtained in §4, that has not yet had the effects of the foreground 
mask accounted for. 

We seek a spectrum Ci such that 

-2 ln£(a,|Cf) = - Cf h ) Q W {C v C^) = -2 ln£(C?|Cf). (C3) 

xv 

where Qw is the inverse-covariance matrix of the combined spectrum which comes with the estimate of Ci . 
Suppose, for simplicity, that £f u n is diagonal, (Sf u n)^, — (Sf u ii)) J (5ii' , then it is straightforward to show that 
the deconvolved, weighted-average spectrum 

Ci = gjjgg^l , (C4) 
Qi 

with 

Q« = 5>!(OM. (C5) 
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is the desired spectrum. Substituting equations (C4) and (C5) into equation (C3) produces equation (C2) 
up to a term which has a weak dependence on comological model, which we ignore. This combined spectrum 
is equivalent to the result we would obtain using the "optimal data compression" of Tcgmark ct al. (1997). 

If the inverse covariance matrix is not diagonal in Z, it can be shown that the optimal combined spectrum 
is given by 

p, v C\ ( S fuii)zz' Cf, d w\, 

Cl = — T^cr — ' (C6) 

with 

Qu> = Y, w \Viuxju>4, (C7) 

ij 

and where Cf d is a fiducial cosmological model which we take to be a flat ACDM model with £lb/i 2 = 0.021, 
Vl c h 2 = 0.129, fitot = 1, h = 0.68, n s = 1.2, and r = 0.2. While this model has parameters values that are 
substantially different than the best- fit WMAP model obtained (afterwards) by Spergel et al. (2003), the 
parameter degeneracies are such that Cf d is close to the best-fit model for I > 100 where this estimator is 
actually used (see §5.1). The combined spectrum is optimal if the fiducial model chosen is the correct one; 
otherwise it is still unbiased but slightly sub-optimal (Gupta & Heavens 2002). 



D. CUT-SKY FISHER MATRIX 

The WMAP sky maps have nearly diagonal pixel-pixel noise covariance (Hinshaw et al. 2003a) which 
greatly simplifies the properties of the power spectrum Fisher matrix. In this Appendix we present an 
analytic derivation of the effect of a sky cut and non-uniform pixel weighting. In D.l, we assume that 
we have an optimal estimator of the power spectrum. In the noise dominated limit, we can obtain an 
exact expression, while in the signal dominated limit, we need to approximate the signal correlation matrix 
to obtain an analytic expression. In D.2, we interpolate the Fisher matrix between the signal and noise- 
clominatcd regimes and show that it agrees with numerical estimates. In D.3, we estimate the power spectrum 
covariance matrix from Monte Carlo simulations of the sky and calibrate the interpolation formula. 



D.l. Cut-Sky Fisher Matrix: Analytic Evaluation 

We equate the Fisher matrix of the power spectrum to the curvature of the likelihood function, equa- 
tion (A16), then develop an approximate form for the covariance matrix, C = S + N, where S is the signal 
matrix and N is the noise matrix. We split the noise matrix into two pieces: a weight term and a mask 
term. In the limit that the pixel noise is diagonal, the weight term has the form 

(Nw 1 )^ = —^Sij = WiSij (Dl) 

CTq 

where i and j are pixel indices, rij is the number of observations of pixel i, ctq is the rms noise of a single 
observation, and, by definition, Wi is the weight of pixel i. The mask, Mi, is defined so that Mi equals 
in pixels that are not used due to foreground contamination, and equals 1 otherwise. The noise matrix can 
then be written as the product of the two terms 

N- 1 = N- X M = MN" 1 = ^8ij = WiSij (D2) 
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where hi = rii and Wi — u>i in the unmasked pixels, and are set to otherwise. We thus define the covariance 
matrix over the full sky, which allows us to exploit the orthogonality properties of the spherical harmonics. 
Note that M 2 = M. 

The inverse of the full covariance matrix can now be written as 

C- 1 = (S + N)- 1 = iNT^SiNT 1 + (D3) 
= M- 1 (SM- 1 +N W )- 1 . (D4) 

The covariance matrix has two limits. In the noise dominated limit, SN _1 <C I, 

cr^isr 1 . (D5) 

In the signal dominated limit, SM _1 <C N w , only the mask alters the covariance matrix, so 

C-^MS^M, (D6) 
where we have set the inverse of the mask to zero where there is no data, i.e. M _1 = M. 



D.l.l. Cut-Sky Fisher Matrix in the Noise Dominated Limit 
In the noise dominated limit, the Fisher matrix, equation (A16), is approximately 

F$ = ha (C-^iC^Pv) - itr (N^N" 1 ^,) . (D7) 

Using 

(Pl)« = ^ = { ^tr- W > P ^ COS %) = W * E Y l*m(Pi) Y lm(Pj) (D8) 

171 

where Oij is the angle between pixels i and j, Y/ m (pi) is a spherical harmonic evaluated in the direction of 
pixel i, and we have employed the addition theorem for spherical harmonics in the last equality. Then the 
Fisher matrix takes the form 



F$ = £ W E E mYCMYimipj) wjY^, {pi)Yv m , (pi). (D9) 
Now expand the weight array as 



2 

mm' ij 



Wi 



--^2wi m Yi m (pi) (D10) 

Ira 

and substitute this into equation (D9) to yield 

F U> = \wiW V E E E ™l"rn"Wl"'m'"^Y l n m n{jp i )Y? m (piW 



mm' l"m" l"'m"' 



(Dll) 

Since Wi was defined over the full sky, the sum over pixels may be expressed in terms of products of Wigner 
3-j symbols. As shown in Appendix E, we may use an orthogonality property of the 3-j symbols to reduce 
the expression for the Fisher matrix to 

Fu ' = 2 WlWl ' % W, "{o J (D12) 

where Wj = J2 m |w; m | 2 , and fl p is the solid angle per pixel. 



D.1.2. Cut-Sky Fisher Matrix in the Signal Dominated Limit 
In the signal dominated limit, the Fisher matrix, equation (A16), is approximately 

Ff v = itr (Cr^CT 1 ^) - itr (s^S" 1 ^) , (D13) 

where S -1 = MiS^ Mj. Using equation (D8) we can write Fp t , in the form 

Ftv = E E MiMjS^j 1 Y* m {Pj)Ylm{Pk) E M k M p S£ Y l*m>(Pp)Yl>m>(Pi). (D14) 

mm' ijk kp 

Substituting a normalized expression for S _1 in the pixel basis, 

S J ^pEc^ YC m { Pi )Y lm { Pi ), (D15) 

lm 

into equation (D14), one obtains a term 

l" l " Wl " m" i j 

Since our mask cuts only ~15% of the sky, the coupling sum, J2i M^Yy m * (pi)Yi?, m ,/(pi) , peaks very sharply 
at \V — l"\ <C Therefore, one may approximate equation (D16) with 

^ ' l"m" i j 

(CiwiGi'Wi' yi A t-r 1 

where, in the last equality, we have used the completeness relation for the spherical harmonics 

E YiMYiM = n; 1 ^-. (D18) 

lm 

With this approximation equation (D14) reduces to 



1 E E^^^^fe^i^fe^'m'W. (D19) 



2 CjCi 

mm' 



Following the same steps outlined above for the derivation of equation (D12) yields 

s _l 1 (2Z + 1)(2Z' + 1) v (I V l"\ 2 



2 CiCv 4tt tjf V 

where M; = J2 m \ m im\ 2 , and m; m is the spherical harmonic transform of the mask. Note that M = 47r/ s 2 ky , 
where / s k y is the fraction of the sky that survives the mask. 
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D.2. Interpolating the Cut-Sky Fisher Matrix 

We can combine the two limiting cases of the Fisher matrix, obtained in the previous section, into a 
single expression 

Fu ' - 2(C l + N l )(C l >+N l ,) Fll '> (D21) 
where TV; = VL p (Jq / '(n bs w i) is the deconvolved noise power spectrum. Using equations (D12) and (D20) 
derived in the previous section, Fu> can be expressed in the low I (signal-dominated) and high I (noise- 
dominated) limits as 

V(2Z + 1)(2Z' + 1) v / I V I" \ 2 

4^ \? Ml "\* Oi' (D22) 



and 



V(2Z + l)(2i' + l)^ Tl . / I V I" 



V Vo 

By construction, these matrices are normalized to 1 on the diagonal, Ff t = F^ = 1, since M = 47r/ s 2 ky and 

We have computed i 7 ';™ ask directly, for selected / values, by evaluating Yi m C~ 1 Yi> m i using the pre- 
conditioner code described by Oh et al. (1999). We find that deviations between the numerical and analytical 
results are consistent with numerical noise in the Fisher matrix estimate. Figure 14 shows the estimated 
Fisher matrix for the Kp2 cut in the noise and the signal dominated limits. 

We interpolate between the two regimes with the following expression: 

Flv ~ ic^mTW) • (D24) 

Note that the form of the Fisher matrix primarily depends on Al, but is weakly dependent on I - there is 
more coupling in the noise dominated limit. 



D.3. Power Spectrum Covariance: Monte Carlo Evaluation 

As discussed in Appendix A. 1.2 we compute the Ci using three different pixel weightings. The uniform 
weighting and the inverse-noise inverse weighting are optimal in the signal-dominated regime and the noise- 
dominated regime, respectively. In between these limits the transitional weighting performs better. In order 
to determine which ranges in I correspond to which regimes, and to calibrate our ansatz for the covariance 
matrix, equation (D21), we proceed as follows. Using 100,000 Monte Carlo simulations of signal plus noise 
(with WMAP noise levels and symmetrized beams), we compute the diagonal elements of the covariance 
matrix for the three different weighting schemes, evaluated with the Kp2 sky cut. Denote these estimates 
Df m 

We find that the uniform weighting produces the smallest D* lm below I = 200. Inverse-noise weighting 
is the best above I = 450, and the transitional weighting produces the lowest variance in between. In each of 
these regimes we use the resulting D; lm to calibrate our ansatz for the diagonal elements of the covariance 
matrix 

D -wrnk lc,+Nt "^ (D25) 
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as illustrated in Verde et al. (2003) . These calibrations produce a smooth correction to equation (D24) of at 
most 6%. No correction at all is required in the signal dominated regime. 



E. SOME USEFUL PROPERTIES OF SPHERICAL HARMONICS 

The derivation of the form of the Fisher matrix in the signal and noise dominated limits led to expressions 
which included a term 

E Y l"™- (PM™ (Pi) y '"'™"' (Pi) Y v m > (Pj)Yl'm' (Pi), (El) 

mm' ij 

where the sum is effectively a double integral over the full sky. This can be evaluated in terms of the Wigner 
3-j symbols, defined as 



dQ n Y im (n)y; ro ,(n)Y^ ro ,,(n) = (-l) m 
Substituting this into equation (El) gives 



(2/ + l)(2/' + l)(2/" + l) / I V V I" 



4tt VOOOMto -ml m" 



(E2) 



J2 (-l) m+m ' {2l + l) ^ V + l) ,J{2l" + l){2V" + l)( 



I V I" \ / I I' V" 

oo o/loo 



\ m —m ml ) \ m — ml ml ) 
We simplify this using the Wigner 3-j orthogonality condition 

y2(-ir +m '( l v , v , l l) = ^-Ts ri , // 8 m ,, m ,,,s(ij , ,n, (E4) 

\ — m ml ml ) \ m — ml ml ) 21" + 1 



where 5(1, V, I") = 1 for |2 — l'\ < I" < I + 1' and is otherwise. This reduces equation (El) to 

(2i + l)(2i' + l) / I V I" s 1 



Att^I \ 



(E5) 



where the factor of fl^ accounts for the fact that equation (El) is a double sum over pixels, instead of a 
double integral. 
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Fig. 1. — The effective noise as a function of I for the 8 differencing assemblies used in the combined power 
spectrum analysis. These spectra were computed from end-end simulations of noise maps, as discussed in 
§2.2. For illustration, the spectra shown here were computed using the quadratic estimator with uniform 
pixel weights. The actual noise model uses three separate weighting schemes in three separate I ranges, and 
thus has discontinuities where the effective noise level changes. The weights are defined in Appendix A. 1.2. 
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Fig. 2. — The full set of individual cross-power spectra for I < 500, computed from the 8 high frequency 
differencing assemblies Ql through W4, 28 spectra in all. The spectra were evaluated from the uncorrected 
sky maps (no Galaxy model subtracted) using the Kp2 sky cut with uniform weighting. The data are plotted 
in color by effective frequency ^JT^Vj with red corresponding to 41 GHz and blue to 94 GHz. The top panel 
shows a very robust measurement of the first acoustic peak with a maximum near / ~ 220. There is also a 
clear indication of the rise to a second peak at I ~ 540 as discussed in §6. The bottom panel shows the ratio of 
each channel to the combined spectrum presented in §5. This clearly shows the residual foreground emission 
due to diffuse Galactic radio emission at low I and to point sources at higher I. The level of contamination, 
which is strongest at Q band, is consistent with the expected level of foreground emission. See Figure 3 for 
the spectra after foreground subtraction. 
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Fig. 3. — The same set of cross-power spectra as shown in Figure 2. Here, the foreground model discussed in 
§3.1 has been subtracted from each channel. The bottom panel shows the ratio of each of the 28 cross-power 
spectra to the combined spectrum presented in §5. Aside from a ~10% discrepancy in the Q band data at 
I < 20, the data are consistent with each other to the sensitivity limits of the individual spectra. Because 
the WMAP data are not sensitivity limited at low I, we use only V and W band data in the final combined 
spectrum for I < 100 (see §5) to minimize residual Galactic contamination. 
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Fig. 4. — The same set of cross-power spectra as shown in Figure 3, but showing the low I spectrum unbinncd. 
The agreement between the individual spectra is striking. The low value of the quadrupole moment, C2, that 
was first seen by COBE-DMR is also seen in the WMAP data. The steep, nearly linear rise in the spectrum 
from I = 2 to 5 translates to a near absence of power in the angular correlation function at separations larger 
than -60° (Spergel et al. 2003; Bennett et al. 2003b). This was also seen in the COBE-DMR data, but it 
is now clear that this is not due to Galaxy modeling errors. The bottom panel shows the fractional rms 
among the 28 WMAP cross-power spectra in black, while the red curve shows the same statistic averaged 
over an ensemble of 1000 Monte Carlo realizations. Based on this we estimate the measurement error on the 
combined spectrum to be <2-5% for I < 100. 
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Fig. 5. — The auto-power spectrum of the combined Q+V+W map evaluated with the three weighting 
schemes defined in Appendix A. 1.2. In each case, the spectrum was computed over the entire I range, and 
black shows uniform weights, red shows inverse- noise weights, and green shows transitional weights. The 
agreement is excellent. 
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Fig. 6. — This figure compares the auto-power spectrum computed from the combined Q+V+W map (grey) 
to the optimally combined cross-power spectrum (black). The close agreement indicates that the noise 
properties of the first-year WMAP data are well understood 
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Fig. 7. — A comparison of the power spectrum computed with data from the ecliptic plane (black) vs. data 
from the ecliptic poles (grey). Note that some of the "bite" features that appear in the combined spectrum 
are not robust to data excision. There is also no evidence that beam ellipticity, which would be more manifest 
in the plane than in the poles, systematically biases the spectrum. This is consistent with estimates of the 
effect given by Page et al. (2003a). 
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Fig. 8. — The final angular power spectrum, 1(1 + l)Ci/2n, obtained from the 28 cross-power spectra, 
as described in §5. The data are plotted with ler measurement errors only which reflect the combined 
uncertainty due to noise, beam, calibration, and source subtraction uncertainties. The solid line shows the 
best-fit ACDM model from Spergel et al. (2003). The grey band around the model is the la uncertainty 
due to cosmic variance on the cut sky. For this plot, both the model and the error band have been binned 
with the same boundaries as the data, but they have been plotted as a splined curve to guide the eye. On 
the scale of this plot the unbinned model curve would be virtually indistinguishable from the binned curve 
except in the vicinity of the third peak. 
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Fig. 9. — A compilation of recent CMB power spectrum measurements compared to the best-fit ACDM 
model from the first-year WMAP data. The data points include noise and cosmic variance uncertainty (but 
not calibration uncertainty) thus we omit the cosmic variance band from the model curve in the Figure. On 
average, the pre- WMAP data agree well with the WMAP power spectrum. 
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Fig. 10. — The WMAP combined power spectrum, in black, compared to a compilation of all CMB data 
published prior to WMAP from Wang et al. (2002a), in red. The WMAP data are plotted with cosmic 
variance plus measurement uncertainties here in order to facilitate a comparison with the compiled data 
which is reported in this way. The data agree well on COBE scales, I < 20, (but note that the WMAP 
cosmic variance errors are computed from the best-fit model rather than the data, thus they appear larger 
than the COBE errors at the quadrupole). However, the overall normalization of the WMAP spectrum is 
~10% higher on smaller scales. 
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Fig. 11. — The WMAP combined power spectrum compared to the locus of predicted spectra, in red, based 
on a joint analysis of pre- WMAP CMB data and 2dFGRS large-scale structure data (Percival et al. 2002). 
As in Figure 8, the WMAP data are plotted with measurement uncertainties, and the best-fit ACDM model 
(Spergel et al. 2003) is plotted with a la cosmic variance error band. The locus of predicted spectra lie 
systematically below the WMAP data at intermediate I. 
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Fig. 12. — A comparison of two power spectra computed from a single V band map. The black points result 
from the maximum likelihood method, the red points from the quadratic estimator computed with uniform 
pixel weighting. In both cases a noise model has been assumed to treat the noise bias. At higher /, the two 
spectra would be nearly identical since both impose inverse noise weighting on the data. 
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Fig. 13. — The quadratic estimator discussed in §A.l has been extensively tested with Monte Carlo sim- 
ulations of first-year WMAP data. The top panel shows a model spectrum in red, and the mean of 600 
realizations of the cross-power spectrum, computed from the V and W band-averaged maps, in black. The 
red error bars show the uncertainty in the mean in bins with Al = 20. The cross-power spectrum estimator 
is unbiased. The bottom panel is similar, but with auto-power spectra computed from the W band-averaged 
map. Here the noise bias term was estimated directly from the high I tail of the computed spectrum. This 
simple auto-power spectrum estimator is unbiased if the noise is white, and is suitable for V band data, but 
see Figure 1 for examples of W band noise spectra. 
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Fig. 14. — Slices of the Fisher (or curvature) matrix normalized as Fw / V Fu Fj} Ji , plotted vs. = I — /', 
from the Fisher matrix for the combined spectrum. Black is I = 30, red is I = 300, green is / = 600. 
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Table 1. WMAP Power Spectrum Noise Modcl a 



Parameter 


Ql 


Q2 


VI 


V2 


Wl 


W2 


W3 


W4 




I < 200 


co(xl0 2 ) 


-6.604 


-6.971 


-6.823 


-6.674 


-7.174 


-7.439 


-7.044 


-7.989 


C!(X10 3 ) 


1.934 


4.244 


2.497 


1.395 


3.060 


3.721 


1.459 


4.820 


c 2 (xl0 4 ) 


-3.798 


-8.521 


-5.203 


-1.483 


-6.021 


-6.901 


-2.677 


-5.870 


c 3 (xl0 5 ) 


2.549 


5.921 


3.715 


0.242 


4.147 


4.540 


1.750 


2.167 


200 <l < 450 


c„(xl0 2 ) 


-6.143 


-6.268 


-6.320 


-6.192 


-6.618 


-6.833 


-6.768 


-6.849 


Cl (xl0 4 ) 


-1.334 


1.341 


-0.800 


-1.351 


0.769 


2.582 


0.885 


2.443 


I > 450 


c„(xl0 2 ) 


-6.202 


-6.167 


-6.338 


-6.264 


-6.564 


-6.704 


-6.688 


-6.772 


Cl (xl0 5 ) 


3.288 


4.085 


3.327 


5.979 


7.342 


12.970 


8.851 


24.430 



a Best-fit coefficients for the noise model in equation (8). The units of the output noise are 
mK 2 , thermodynamic temperature. 



